In [1]:
%pylab inline
import matplotlib.pyplot as plt
import collections
from IPython.core.display import HTML
In [2]:
HTML('''<script>
code_show=true;
function code_toggle() {
if (code_show){
$('div.input').hide();
} else {
$('div.input').show();
}
code_show = !code_show
}
$( document ).ready(code_toggle);
</script>
<form action="javascript:code_toggle()"><input type="submit"
value="Click here to toggle on/off the raw code."></form>''')
Out[2]:
In [3]:
# Az abra kimentesehez az alabbiakat a plt.show() ele kell tenni!!!
#savefig('fig_rainbow_p2_1ray.pdf'); # Abra kimentese
#savefig('fig_rainbow_p2_1ray.eps'); # Abra kimentese
# Abra es fontmeretek
xfig_meret= 8 # 12 nagy abrahoz
yfig_meret= 6 # 12 nagy abrahoz
xyticks_meret= 15 # 20 nagy abrahoz
xylabel_meret= 20 # 30 nagy abrahoz
legend_meret= 21 # 30 nagy abrahoz
In [4]:
def free_en(kp, Gvec):
# Gvec ---> G vectors
eig_en = []
for G in Gvec:
eig_en.append(dot((kp-G),(kp-G)))
return sort(eig_en)
In [5]:
def Ham_op(kv, Gvec, params):
# Dirac-delta szennyezok a racsatomokon
# Hsize = the size of the ham matrix, input
w = params # w ----> a Dirac-delta potencial erossege, input
Hsize = len(Gvec)
Ham = w*ones((Hsize,Hsize)) # a Ham Hamiltonian matrix minden eleme = w
# a Ham diagonalis elemei = a kinetic energy:
i=0
for g in Gvec:
Ham[i,i] = dot(kv-g,kv-g)
i +=1
return Ham
Simple cubic lattice¶
In [6]:
## simple cubic lattice
Npoints = 50; ## hany pontbol keszul a gorbe
ymax = 3.1;
## unit cell vectors of the reciprocal vectors:
## in units of 2 pi/a
b1 = array([1,0,0])
b2 = array([0,1,0])
b3 = array([0,0,1])
## high symmetry points in the Brillouin zone
## FONTOS: a fenti abra szerinti pontokat kell venni,
## a szomszedos spec. k pontokat kell venni.
Gp = array([0,0,0])
Xp = 0.5*array([0,1,0])
Rp = 0.5*array([1,1,1])
Mp = 0.5*array([1,1,0])
# creating the reciprocal vectors:
Gnum = 2 #### (2*Gnum+1)**3, number of Gvec
indx=array(range(-Gnum,Gnum+1))
Gvec = []
for n1 in indx:
for n2 in indx:
for n3 in indx:
Gvec.append(n1*b1+n2*b2+n3*b3)
specNev_map = {tuple(Gp):"$\Gamma$",tuple(Xp):"$X$",tuple(Mp):"$M$",
tuple(Rp):"$R$"}
figsize(xfig_meret,yfig_meret)
## G-X-M-G-R-X line
specpoints = array([Gp,Xp,Mp,Gp,Rp,Xp])
klist = []
specpos = [0]
kk = [0]
tmp = 0
for i in (range(len(specpoints)-1)):
sl = specpoints[i+1]-specpoints[i]
sl_hossz = sqrt(dot(sl,sl))
dk = sl/Npoints
for num in arange(0,Npoints):
# k vector along the line given by speclines:
klist.append(specpoints[i] + num*dk)
# length of the k vector and shifted:
tmp += sqrt(dot(dk,dk))
kk.append(tmp)
specpos.append(tmp)
klist.append(klist[-1]+dk)
enlist = []
for i in klist:
enlist.append(free_en(i, Gvec))
enlist= array(enlist)
eigszam= len(enlist[0,:])
print("The # of eigenvalues = ",eigszam,"= # of G vectors =",len(Gvec))
#specNev = ["$\Gamma$","$X$","$M$","$\Gamma$","$R$","$X$"]
# drawing the figure
figsize(xfig_meret,yfig_meret)
for ii in range(0,len(enlist[0,:])):
plot(kk,enlist[:,ii],color='r')
i = 0
for xc in specpos:
axvline(x=xc,color='b')
text(xc*0.95,-0.2,specNev_map[tuple(specpoints[i])],fontsize=xylabel_meret)
i=i+1
### Brillouin zone Volume:
VBZ= dot(b3,cross(b1,b2))
print("Volume of the Brilouin zone = ", VBZ,"\n")
### Fermi energy:
EF_vec = []
for n in range(1,9):
tmp = (3/8/pi*VBZ*n)**(2/3)
EF_vec.append(tmp)
print("Electrons per unit cell = ", n, ", Fermi energy = ",round(tmp,3))
## horizontal line at Fermi energies:
i=1
for ef in EF_vec:
plt.axhline(y=ef, color='b',ls='--')
nel = str(i)
text(specpos[-1]*1.05,ef*0.97,nel,fontsize=10)
i +=1
ylabel(r'$E(k)$',fontsize=xylabel_meret)
# kikapcsoljuk a xticks:
xticks([], [])
xlim(0,specpos[-1])
ylim(0,ymax)
title("simple cubic",fontsize=xylabel_meret)
grid();
#savefig('fifi.eps'); # Abra kimentese
In [7]:
## calculation of degeneracy
print("\nDegeneracy at the special k points without Dirac delta potential: \n=================================\n")
i=0
for sp in specpoints:
se = free_en(sp, Gvec)
sajatertekek =[ round(x,2) for x in se]
deg=collections.Counter(sajatertekek);
keys = sorted(deg.keys())
print("\nDegeneracies in the %s point:" % specNev_map[tuple(specpoints[i])])
for element in keys:
print("Energy: %.2f, degeneracy: %d" % ( element, deg[element]))
i +=1
print("\n")
Simple cubic lattice with Dirac delta potential¶
In [8]:
## simple cubic lattice
Npoints = 50; ## hany pontbol keszul a gorbe
ymax = 2.1;
w = -0.01 #### w =-0.1 ---> a Dirac-delta potencial erossege, input
params=[w]
print("The # of G vectors = (2*Gnum+1)**3 = ", (2*Gnum+1)**3)
print("w = ",w)
## unit cell vectors of the reciprocal vectors:
## in units of 2 pi/a
b1 = array([1,0,0])
b2 = array([0,1,0])
b3 = array([0,0,1])
## high symmetry points in the Brillouin zone
## FONTOS: a fenti abra szerinti pontokat kell venni,
## a szomszedos spec. k pontokat kell venni.
Gp = array([0,0,0])
Xp = 0.5*array([0,1,0])
Rp = 0.5*array([1,1,1])
Mp = 0.5*array([1,1,0])
# creating the reciprocal vectors:
Gnum = 1 #### (2*Gnum+1)**3, number of Gvec
indx=array(range(-Gnum,Gnum+1))
Gvec = []
for n1 in indx:
for n2 in indx:
for n3 in indx:
Gvec.append(n1*b1+n2*b2+n3*b3)
specNev_map = {tuple(Gp):"$\Gamma$",tuple(Xp):"$X$",tuple(Mp):"$M$",
tuple(Rp):"$R$"}
figsize(xfig_meret,yfig_meret)
## G-X-M-G-R-X line
specpoints = array([Gp,Xp,Mp,Gp,Rp,Xp])
klist = []
specpos = [0]
kk = [0]
tmp = 0
for i in (range(len(specpoints)-1)):
sl = specpoints[i+1]-specpoints[i]
sl_hossz = sqrt(dot(sl,sl))
dk = sl/Npoints
for num in arange(0,Npoints):
# k vector along the line given by speclines:
klist.append(specpoints[i] + num*dk)
# length of the k vector and shifted:
tmp += sqrt(dot(dk,dk))
kk.append(tmp)
specpos.append(tmp)
klist.append(klist[-1]+dk)
enlist = []
for kv in klist:
enlist.append(free_en(kv, Gvec))
enlist= array(enlist)
enlistw = []
for kv in klist:
enlistw.append(eigvalsh(Ham_op(kv, Gvec, params)))
enlistw= array(enlistw)
eigszam= len(enlist[0,:])
print("The # of eigenvalues = ",eigszam,"= # of G vectors =",len(Gvec))
specNev = ["$\Gamma$","$X$","$M$","$\Gamma$","$R$","$X$"]
# drawing the figure
figsize(10,5)
subplot(121)
for ii in range(0,len(enlist[0,:])):
plot(kk,enlist[:,ii],color='r')
i = 0
for xc in specpos:
axvline(x=xc,color='b')
text(xc,-0.2,specNev_map[tuple(specpoints[i])],fontsize=xylabel_meret)
i=i+1
ylabel(r'$E(k)$',fontsize=xylabel_meret)
# kikapcsoljuk a xticks:
xticks([], [])
xlim(0,specpos[-1])
ylim(0,ymax)
cim = "sc, w="+str(0)
title(cim,fontsize=xylabel_meret)
grid();
subplot(122)
for ii in range(0,len(enlist[0,:])):
plot(kk,enlistw[:,ii],color='r')
i = 0
for xc in specpos:
axvline(x=xc,color='b')
text(xc,-0.2,specNev[i],fontsize=xylabel_meret)
i=i+1
#ylabel(r'$E(k)$',fontsize=xylabel_meret)
# kikapcsoljuk a xticks:
xticks([], [])
xlim(0,specpos[-1])
ylim(-0.0,ymax)
cim = "sc, w="+str(w)
title(cim,fontsize=xylabel_meret)
grid();
In [9]:
#sort(free_en(Gp, Gvec)), eigvalsh(Ham_op(Gp, Gvec, params))
Special points in the Brillouin zone (in units of $2\pi/a$):¶
| Point | Cartesian coordinates (in units of $2\pi/a$) |
|---|---|
| G | array([0,0,0]) |
| X | array([0,1,0]) |
| W | array([1/2,1,0]) |
| K | array([3/4,3/4,0]) |
| L | array([1/2,1/2,1/2]) |
| U | array([1/4,1,1/4]) |
U,K points together !!!!!!!!¶
Sólyom Jenő: A modern szilárdtest-fizika alapjai II.
Fémek, félvezetők, szupravezetők, 881. old., 20.1 ábra
In [10]:
## fcc lattice:
Npoints = 70; ## hany pontbol keszul a gorbe
ymax = 4.5;
## unit cell vectors of the reciprocal vectors:
## in units of 2 pi/a
b1 = array([-1,1,1])
b2 = array([1,-1,1])
b3 = array([1,1,-1])
## high symmetry points in the Brillouin zone
## FONTOS: a szomszedos spec. k pontokat kell venni
## (ezek most nem a fenti abranak megfelelo pontok, de szomszedos pontok,
## es igy jo az abra)
Gp = array([0,0,0])
Xp = array([0,1,0])
Wp = array([1/2,1,0])
Kp = array([3/4,3/4,0])
Lp = array([1/2,1/2,1/2])
Up = array([1/4,1,1/4])
# creating the reciprocal vectors:
Gnum = 2 #### (2*Gnum+1)**3, number of Gvec
indx=array(range(-Gnum,Gnum+1))
Gvec = []
for n1 in indx:
for n2 in indx:
for n3 in indx:
Gvec.append(n1*b1+n2*b2+n3*b3)
#figsize(xfig_meret,yfig_meret)
figsize(6,6)
## G-X-W-L-G-K-U-X line
#specpoints = array([Gp,Xp,Wp,Lp,Gp,Kp,Up,Xp])
#specpoints = array([Lp,Gp,Xp,Up,Kp,Gp,Wp])
specpoints = array([Lp,Gp,Xp,Up,Kp,Gp]) # Solyom konyvben: II. kotet, 20.1. abra
# Az U,K pontok egybe!!!!!
## FIGYELEM: fcc racs eseten az irodalomban az U,K jeloles azt jelenti, hogy
## U-K vonal menten NEM rajzoljuk ki a savot!!!!!
### kiiratashoz, az x-tengelyen:
specNev_map = {tuple(Gp):"$\Gamma$",tuple(Xp):"$X$",tuple(Wp):"$W$",
tuple(Lp):"$L$",tuple(Up):"$U,K$",tuple(Kp):"$\Gamma$"}
klist = []
specpos = [0]
kk = [0]
tmp = 0
for i in (range(len(specpoints)-1)):
if not i==3: ########## Az U,K pontok egybe!!!!!
sl = specpoints[i+1]-specpoints[i]
sl_hossz = sqrt(dot(sl,sl))
dk = sl/Npoints
for num in arange(0,Npoints):
# k vector along the line given by speclines:
klist.append(specpoints[i] + num*dk)
# length of the k vector and shifted:
tmp += sqrt(dot(dk,dk))
kk.append(tmp)
specpos.append(tmp)
klist.append(klist[-1]+dk)
enlist = []
for i in klist:
enlist.append(free_en(i, Gvec))
enlist= array(enlist)
eigszam= len(enlist[0,:])
print("The # of eigenvalues = ",eigszam,"= # of G vectors =",len(Gvec))
# drawing the figure
figsize(xfig_meret,yfig_meret)
for ii in range(0,len(enlist[0,:])):
plot(kk,enlist[:,ii],color='r')
i = 0
for xc in specpos:
plt.axvline(x=xc,color='b')
text(xc*0.95,-0.35,specNev_map[tuple(specpoints[i])],fontsize=xylabel_meret)
i=i+1
### Brillouin zone Volume:
VBZ= dot(b3,cross(b1,b2))
print("Volume of the Brilouin zone = ", VBZ,"\n")
### Fermi energy:
EF_vec = []
for n in range(1,9):
tmp = (3/8/pi*VBZ*n)**(2/3)
EF_vec.append(tmp)
print("Electrons per unit cell = ", n, ", Fermi energy = ",round(tmp,3))
## horizontal line at Fermi energies:
i=1
for ef in EF_vec:
plt.axhline(y=ef, color='b',ls='--')
nel = str(i)
text(specpos[-1]*1.05,ef*0.97,nel,fontsize=10)
i +=1
ylabel(r'$E(k)$',fontsize=xylabel_meret)
# kikapcsoljuk a xticks:
xticks([], [])
xlim(0,specpos[-1])
ylim(0,ymax)
title("fcc cubic",fontsize=xylabel_meret);
grid();
#savefig('Fig_free_LGXU-KG.png'); # Abra kimentese
In [11]:
## calculation of degeneracy
print("\nDegeneracy at the special k points without Dirac delta potential: \n=================================\n")
i=0
for sp in specpoints:
se = free_en(sp, Gvec)
sajatertekek =[ round(x,2) for x in se]
deg=collections.Counter(sajatertekek);
keys = sorted(deg.keys())
print("\nDegeneracies in the %s point:" % specNev_map[tuple(specpoints[i])])
for element in keys:
print("Energy: %.2f, degeneracy: %d" % ( element, deg[element]))
i +=1
print("\n")
Results below are the same as in Ashcroft & Mermin: Solid State Physics¶
https://www.amazon.com/Solid-State-Physics-Neil-Ashcroft/dp/0030839939
In [12]:
## fcc lattice:
# Ashcroft & Mermin: Solid State Physics
Npoints = 70; ## hany pontbol keszul a gorbe
ymax = 4.5;
## unit cell vectors of the reciprocal vectors:
## in units of 2 pi/a
b1 = array([-1,1,1])
b2 = array([1,-1,1])
b3 = array([1,1,-1])
## high symmetry points in the Brillouin zone
## FONTOS: a szomszedos spec. k pontokat kell venni
## (ezek most nem a fenti abranak megfelelo pontok, de szomszedos pontok,
## es igy jo az abra)
Gp = array([0,0,0])
Xp = array([0,1,0])
Wp = array([1/2,1,0])
Kp = array([3/4,3/4,0])
Lp = array([1/2,1/2,1/2])
Up = array([1/4,1,1/4])
# creating the reciprocal vectors:
Gnum = 2 #### (2*Gnum+1)**3, number of Gvec
indx=array(range(-Gnum,Gnum+1))
Gvec = []
for n1 in indx:
for n2 in indx:
for n3 in indx:
Gvec.append(n1*b1+n2*b2+n3*b3)
#figsize(xfig_meret,yfig_meret)
figsize(6,6)
## G-X-W-L-G-K-U-X line
specpoints = array([Gp,Xp,Wp,Lp,Gp,Kp,Up,Xp])
#specpoints = array([Lp,Gp,Xp,Up,Kp,Gp,Wp])
#specpoints = array([Lp,Gp,Xp,Up,Kp,Gp]) # Solyom konyvben: II. kotet, 20.1. abra
# Az U,K pontok egybe!!!!!
## FIGYELEM: fcc racs eseten az irodalomban az U,K jeloles azt jelenti, hogy
## U-K vonal menten NEM rajzoljuk ki a savot!!!!!
### kiiratashoz, az x-tengelyen:
#specNev_map = {tuple(Gp):"$\Gamma$",tuple(Xp):"$X$",tuple(Wp):"$W$",tuple(Lp):"$L$",tuple(Up):"$U,K$",tuple(Kp):"$\Gamma$"}
specNev_map = {tuple(Gp):"$\Gamma$",tuple(Xp):"$X$",tuple(Wp):"$W$",
tuple(Lp):"$L$",tuple(Gp):"$\Gamma$",
tuple(Kp):"$K$",tuple(Up):"$X$"}
klist = []
specpos = [0]
kk = [0]
tmp = 0
for i in (range(len(specpoints)-1)):
if not i==5: ########## Az U,K pontok egybe!!!!!
sl = specpoints[i+1]-specpoints[i]
sl_hossz = sqrt(dot(sl,sl))
dk = sl/Npoints
for num in arange(0,Npoints):
# k vector along the line given by speclines:
klist.append(specpoints[i] + num*dk)
# length of the k vector and shifted:
tmp += sqrt(dot(dk,dk))
kk.append(tmp)
specpos.append(tmp)
klist.append(klist[-1]+dk)
enlist = []
for i in klist:
enlist.append(free_en(i, Gvec))
enlist= array(enlist)
eigszam= len(enlist[0,:])
print("The # of eigenvalues = ",eigszam,"= # of G vectors =",len(Gvec))
## calculation of degeneracy
# drawing the figure
figsize(xfig_meret,yfig_meret)
for ii in range(0,len(enlist[0,:])):
plot(kk,enlist[:,ii],color='r')
i = 0
for xc in specpos:
plt.axvline(x=xc,color='b')
text(xc,-0.35,specNev_map[tuple(specpoints[i])],fontsize=xylabel_meret)
i=i+1
### Brillouin zone Volume:
VBZ= dot(b3,cross(b1,b2))
print("Volume of the Brilouin zone = ", VBZ,"\n")
### Fermi energy:
EF_vec = []
for n in range(1,9):
tmp = (3/8/pi*VBZ*n)**(2/3)
EF_vec.append(tmp)
print("Electrons per unit cell = ", n, ", Fermi energy = ",round(tmp,3))
## horizontal line at Fermi energies:
i=1
for ef in EF_vec:
plt.axhline(y=ef, color='b',ls='--')
nel = str(i)
text(specpos[-1]*1.05,ef*0.97,nel,fontsize=10)
i +=1
ylabel(r'$E(k)$',fontsize=xylabel_meret)
# kikapcsoljuk a xticks:
xticks([], [])
xlim(0,specpos[-1])
ylim(0,ymax)
title("fcc, $K,U$ pontok egybe",fontsize=xylabel_meret);
grid();
#savefig('Fig_free_GXWLGK-UX.png'); # Abra kimentese
In [13]:
## simple fcc lattice
Npoints = 50; ## hany pontbol keszul a gorbe
ymax = 4.2;
w = -0.02 #### w ---> a Dirac-delta potencial erossege, input
params=[w]
print("The # of G vectors = (2*Gnum+1)**3 = ", (2*Gnum+1)**3)
print("w = ",w)
## unit cell vectors of the reciprocal vectors:
## in units of 2 pi/a
b1 = array([-1,1,1])
b2 = array([1,-1,1])
b3 = array([1,1,-1])
## high symmetry points in the Brillouin zone
## FONTOS: a szomszedos spec. k pontokat kell venni
## (ezek most nem a fenti abranak megfelelo pontok, de szomszedos pontok,
## es igy jo az abra)
Gp = array([0,0,0])
Xp = array([0,1,0])
Wp = array([1/2,1,0])
Kp = array([3/4,3/4,0])
Lp = array([1/2,1/2,1/2])
Up = array([1/4,1,1/4])
# creating the reciprocal vectors:
Gnum = 2 #### (2*Gnum+1)**3, number of Gvec
indx=array(range(-Gnum,Gnum+1))
Gvec = []
for n1 in indx:
for n2 in indx:
for n3 in indx:
Gvec.append(n1*b1+n2*b2+n3*b3)
specNev_map = {tuple(Gp):"$\Gamma$",tuple(Xp):"$X$",tuple(Wp):"$W$",
tuple(Lp):"$L$",tuple(Kp):"$K$",tuple(Up):"$U$"}
figsize(xfig_meret,yfig_meret)
## G-X-W-L-G-K-X line
#specpoints = array([Gp,Xp,Wp,Lp,Gp,Kp,Xp])
#specpoints = array([Lp,Gp,Xp,Up,Kp,Gp,Wp])
specpoints = array([Lp,Gp,Xp,Up,Kp,Gp]) # Solyom konyvben: II. kotet, 20.1. abra
# Az U,K pontok egybe!!!!!
## FIGYELEM: fcc racs eseten az irodalomban az U,K jeloles azt jelenti, hogy
## U-K vonal menten NEM rajzoljuk ki a savot!!!!!
### kiiratashoz, az x-tengelyen:
specNev_map = {tuple(Gp):"$\Gamma$",tuple(Xp):"$X$",tuple(Wp):"$W$",
tuple(Lp):"$L$",tuple(Up):"$U,K$",tuple(Kp):"$\Gamma$"}
klist = []
specpos = [0]
kk = [0]
tmp = 0
for i in (range(len(specpoints)-1)):
if not i==3: ########## Az U,K pontok egybe!!!!!
sl = specpoints[i+1]-specpoints[i]
sl_hossz = sqrt(dot(sl,sl))
dk = sl/Npoints
for num in arange(0,Npoints):
# k vector along the line given by speclines:
klist.append(specpoints[i] + num*dk)
# length of the k vector and shifted:
tmp += sqrt(dot(dk,dk))
kk.append(tmp)
specpos.append(tmp)
klist.append(klist[-1]+dk)
enlist = []
for kv in klist:
enlist.append(free_en(kv, Gvec))
enlist= array(enlist)
enlistw = []
for kv in klist:
enlistw.append(eigvalsh(Ham_op(kv, Gvec, params)))
enlistw= array(enlistw)
eigszam= len(enlist[0,:])
print("The # of eigenvalues = ",eigszam,"= # of G vectors =",len(Gvec))
# drawing the figure
figsize(10,5)
subplot(121)
for ii in range(0,len(enlist[0,:])):
plot(kk,enlist[:,ii],color='r')
i = 0
for xc in specpos:
axvline(x=xc,color='b')
text(xc-0.1,-0.3,specNev_map[tuple(specpoints[i])],fontsize=xylabel_meret)
i=i+1
ylabel(r'$E(k)$',fontsize=xylabel_meret)
# kikapcsoljuk a xticks:
xticks([], [])
xlim(0,specpos[-1])
ylim(0,ymax)
cim = "fcc, w="+str(0)
title(cim,fontsize=xylabel_meret)
grid();
subplot(122)
for ii in range(0,len(enlist[0,:])):
plot(kk,enlistw[:,ii],color='r')
i = 0
for xc in specpos:
axvline(x=xc,color='b')
text(xc-0.1,-0.3,specNev_map[tuple(specpoints[i])],fontsize=xylabel_meret)
i=i+1
#ylabel(r'$E(k)$',fontsize=xylabel_meret)
# kikapcsoljuk a xticks:
xticks([], [])
xlim(0,specpos[-1])
ylim(0,ymax)
cim = "fcc, w="+str(w)
title(cim,fontsize=xylabel_meret)
grid();
#savefig('Fig_free_quasi_LGXU-KG.png'); # Abra kimentese
In [14]:
free_en([0,0,0], Gvec)
Out[14]:
In [15]:
## Az K es U points nem teljesen azonosak:
free_en(Kp, Gvec),free_en(Up, Gvec)
Out[15]:
In [16]:
## bcc lattice:
Npoints = 50; ## hany pontbol keszul a gorbe
ymax = 1.4;
## unit cell vectors of the reciprocal vectors:
## in units of 2 pi/a
b1 = 0.5*array([0,1,1])
b2 = 0.5*array([1,0,1])
b3 = 0.5*array([1,1,0])
## high symmetry points in the Brillouin zone
## FONTOS: a szomszedos spec. k pontokat kell venni
## (ezek most nem a fenti abranak megfelelo pontok, de szomszedos pontok,
## es igy jo az abra)
Gp = array([0,0,0])
Hp = 0.5*array([0,1,0])
Pp = 0.25*array([1,1,1])
Np = 0.25*array([1,1,0])
# creating the reciprocal vectors:
Gnum = 2 #### (2*Gnum+1)**3, number of Gvec
indx=array(range(-Gnum,Gnum+1))
Gvec = []
for n1 in indx:
for n2 in indx:
for n3 in indx:
Gvec.append(n1*b1+n2*b2+n3*b3)
specNev_map = {tuple(Gp):"$\Gamma$",tuple(Hp):"$H$",tuple(Pp):"$P$",
tuple(Np):"$N$"}
figsize(xfig_meret,yfig_meret)
## G-H-N-G-P-H line
specpoints = array([Gp,Hp,Np,Gp,Pp,Hp])
klist = []
specpos = [0]
kk = [0]
tmp = 0
for i in (range(len(specpoints)-1)):
sl = specpoints[i+1]-specpoints[i]
sl_hossz = sqrt(dot(sl,sl))
dk = sl/Npoints
for num in arange(0,Npoints):
# k vector along the line given by speclines:
klist.append(specpoints[i] + num*dk)
# length of the k vector and shifted:
tmp += sqrt(dot(dk,dk))
kk.append(tmp)
specpos.append(tmp)
klist.append(klist[-1]+dk)
enlist = []
for i in klist:
enlist.append(free_en(i, Gvec))
enlist= array(enlist)
eigszam= len(enlist[0,:])
print("The # of eigenvalues = ",eigszam,"= # of G vectors =",len(Gvec))
# drawing the figure
figsize(xfig_meret,yfig_meret)
for ii in range(0,len(enlist[0,:])):
plot(kk,enlist[:,ii],color='r')
i = 0
for xc in specpos:
plt.axvline(x=xc,color='b')
text(xc*0.95,-0.1,specNev_map[tuple(specpoints[i])],fontsize=xylabel_meret)
i=i+1
### Brillouin zone Volume:
VBZ= dot(b3,cross(b1,b2))
print("Volume of the Brilouin zone = ", VBZ,"\n")
### Fermi energy:
EF_vec = []
for n in range(1,9):
tmp = (3/8/pi*VBZ*n)**(2/3)
EF_vec.append(tmp)
print("Electrons per unit cell = ", n, ", Fermi energy = ",round(tmp,3))
## horizontal line at Fermi energies:
i=1
for ef in EF_vec:
plt.axhline(y=ef, color='b',ls='--')
nel = str(i)
text(specpos[-1]*1.05,ef*0.97,nel,fontsize=10)
i +=1
ylabel(r'$E(k)$',fontsize=xylabel_meret)
# kikapcsoljuk a xticks:
xticks([], [])
xlim(0,specpos[-1])
ylim(0,ymax)
title("bcc cubic",fontsize=xylabel_meret)
grid();
In [17]:
## calculation of degeneracy
print("\nDegeneracy at the special k points without Dirac delta potential: \n=================================\n")
i=0
for sp in specpoints:
se = free_en(sp, Gvec)
sajatertekek =[ round(x,2) for x in se]
deg=collections.Counter(sajatertekek);
keys = sorted(deg.keys())
print("\nDegeneracies in the %s point:" % specNev_map[tuple(specpoints[i])])
for element in keys:
print("Energy: %.2f, degeneracy: %d" % ( element, deg[element]))
i +=1
print("\n")
In [18]:
## simple bcc lattice
Npoints = 50; ## hany pontbol keszul a gorbe
ymax = 0.8;
w = -0.005 #### w ---> a Dirac-delta potencial erossege, input
params=[w]
print("The # of G vectors = (2*Gnum+1)*(2*Gnum+1) = ", (2*Gnum+1)*(2*Gnum+1))
print("w = ",w)
## unit cell vectors of the reciprocal vectors:
## in units of 2 pi/a
b1 = 0.5*array([0,1,1])
b2 = 0.5*array([1,0,1])
b3 = 0.5*array([1,1,0])
## high symmetry points in the Brillouin zone
## FONTOS: a szomszedos spec. k pontokat kell venni
## (ezek most nem a fenti abranak megfelelo pontok, de szomszedos pontok,
## es igy jo az abra)
Gp = array([0,0,0])
Hp = 0.5*array([0,1,0])
Pp = 0.25*array([1,1,1])
Np = 0.25*array([1,1,0])
# creating the reciprocal vectors:
Gnum = 2 #### (2*Gnum+1)**3, number of Gvec
indx=array(range(-Gnum,Gnum+1))
Gvec = []
for n1 in indx:
for n2 in indx:
for n3 in indx:
Gvec.append(n1*b1+n2*b2+n3*b3)
figsize(xfig_meret,yfig_meret)
## G-H-N-G-P-H line
specpoints = array([Gp,Hp,Np,Gp,Pp,Hp])
klist = []
specpos = [0]
kk = [0]
tmp = 0
for i in (range(len(specpoints)-1)):
sl = specpoints[i+1]-specpoints[i]
sl_hossz = sqrt(dot(sl,sl))
dk = sl/Npoints
for num in arange(0,Npoints):
# k vector along the line given by speclines:
klist.append(specpoints[i] + num*dk)
# length of the k vector and shifted:
tmp += sqrt(dot(dk,dk))
kk.append(tmp)
specpos.append(tmp)
klist.append(klist[-1]+dk)
enlist = []
for kv in klist:
enlist.append(free_en(kv, Gvec))
enlist= array(enlist)
enlistw = []
for kv in klist:
enlistw.append(eigvalsh(Ham_op(kv, Gvec, params)))
enlistw= array(enlistw)
eigszam= len(enlist[0,:])
print("The # of eigenvalues = ",eigszam,"= # of G vectors =",len(Gvec))
# drawing the figure
figsize(10,5)
subplot(121)
for ii in range(0,len(enlist[0,:])):
plot(kk,enlist[:,ii],color='r')
i = 0
for xc in specpos:
axvline(x=xc,color='b')
text(xc,-0.05,specNev_map[tuple(specpoints[i])],fontsize=xylabel_meret)
i=i+1
ylabel(r'$E(k)$',fontsize=xylabel_meret)
# kikapcsoljuk a xticks:
xticks([], [])
xlim(0,specpos[-1])
ylim(0,ymax)
cim = "bcc, w="+str(0)
title(cim,fontsize=xylabel_meret)
grid();
subplot(122)
for ii in range(0,len(enlist[0,:])):
plot(kk,enlistw[:,ii],color='r')
i = 0
for xc in specpos:
axvline(x=xc,color='b')
text(xc,-0.05,specNev_map[tuple(specpoints[i])],fontsize=xylabel_meret)
i=i+1
ylabel(r'$E(k)$',fontsize=xylabel_meret)
# kikapcsoljuk a xticks:
xticks([], [])
xlim(0,specpos[-1])
ylim(0,ymax)
cim = "bcc, w="+str(w)
title(cim,fontsize=xylabel_meret)
grid();
In [19]:
#tick_params(
# axis='x', # changes apply to the x-axis
# which='both', # both major and minor ticks are affected
# bottom='off', # ticks along the bottom edge are off
# top='off', # ticks along the top edge are off
# labelbottom='off') # labels along the bottom edge are off
In [ ]:
In [ ]: